R Notebook

This file will take you through recreating the figures presented in the manuscript. Note that some of the figures went through additional edits in Illustrator for additional labelling!

Preparation

Load packages

library(splancs)

Loading required package: sp
Warning: package ‘sp’ was built under R version 4.2.3

Spatial Point Pattern Analysis Code in S-Plus
 
 Version 2 - Spatial and Space-Time analysis


Attaching package: ‘splancs’

The following object is masked from ‘package:tidyr’:

    tribble

The following object is masked from ‘package:dplyr’:

    tribble

load variables

parameters_tsukushi <- c(R1 = 8.89*10^6, 
                lambda = 3.7*10^5,
                mu = 0.025, 
                p = 8*10^-6, 
                alpha = 1, 
                alphag = 2, 
                beta = 5.721, 
                mum = 48, 
                mug = 4, 
                I0 = 43.85965, 
                Ig0 = 0, 
                a = 150, 
                b = 100, 
                sp = 1,
                psin = 16.69234,
                psiw = 0.8431785,
                phin = 0.03520591, 
                phiw = 550.842,
                iota = 2.18*(10^6),
                rho = 0.2627156)

# import in data files
ez_label <- read.csv(here("code_repository/data/ez_label.csv")) ## labelling schene
si_opt.df <- read.csv(here("code_repository/data/si_opt.csv")) ## optimized parameter + fitness list
si_dyn.df <- read_parquet(here("code_repository/data/si_dyn.parquet")) ## dynamics of single cue models
si_rn.df <-  read_parquet(here("code_repository/data/si_rn.parquet")) ## reaction norms of single cue models
si_rug.df <- read_parquet(here("code_repository/data/si_rug.parquet")) ## data for rug plots
mc_all_fitness.df <- read_parquet(here("code_repository/data/mc_all_fitness.parquet")) ## fitness values for mc when all parameters are varying
mc_single_fitness.df <- read_parquet(here("code_repository/data/mc_single_fitness.parquet")) ## fitness values for mc when only one parameter is varying
dual_cue_f_lc.df <- read.csv(here("code_repository/data/dual_cue_fitness_local.csv")) ## optimized dual cue strategy using L-BFGS-B
dual_cue_f_glb.df <- read.csv(here("code_repository/data/dual_cue_fitness_global.csv")) ## optimized dual cue strategy using DEoptim + L-BFG
dual_cue_f_final.df <- read.csv(here("code_repository/data/dual_cue_fitness_final.csv")) ## final dataframe containing strategy that produced highest fitness (dual cue models)
dual_selected_cr.df <- read_parquet(here("code_repository/data/dual_selected_cr.parquet")) ## cr dynamics of selected models
cue_range_si_alt.df <- read.csv(here("code_repository/data/cue_range_si_alt.csv")) ## cue ranges for dual cue models
dual_cue_dyn.df <- read_parquet(here("code_repository/data/dual_cue_dyn.parquet")) ## dynamics of dual cue models
exp_ss.df <- read.csv(here("code_repository/data/experimental_data.csv")) ## cleaned experimental records of P. chabaudi infection
posterior.df <- read.csv(here("code_repository/data/posterior.csv")) ## posterior distribution of different parameter values
mc_burst.df <- read.csv(here("code_repository/data/mc_burst_opt.csv")) ## fitness values of cues when optimized with different beta values
mc_burst_single_input.df <- read.csv(here("code_repository/data/mc_burst_single_input.csv")) ## iput dataframe for simulating dyamics for for models optimized with different burst values (single cues only)
mc_burst_dual_input.df <- read.csv(here("code_repository/data/mc_burst_dual_input.csv")) ## ditto for R log + I log
mc_burst_final_diff.df <- read.csv(here("code_repository/data/mc_burst_final_diff.csv")) ## fitness of parasites adopting different parameters optimized for diff burst size but simulated with discordant burst size 
mc_burst_rn.df <- read_parquet(here("code_repository/data/mc_burst_rn.parquet")) ## dataframe containing reaction norms of strains optimzied at various burst size filtered by relevant ranges (rug)
mc_burst_rn_RI.df <- read_parquet(here("code_repository/data/mc_burst_rn_RI.parquet")) ## dataframe containing reaction norms of strains optimzied at various burst size filtered by relevant ranges (rug)
mc_burst_rug.df_f <- read_parquet(here("code_repository/data/mc_burst_rug_f.parquet")) ## rug dataframe for parasites adopting the ideal strategy for various burst sizes
mc_burst_RI_rug.df <- read_parquet(here("code_repository/data/mc_burst_RI_rug.parquet")) ## dito but for R log and I log
mc_burst_dyn.df <- read_parquet(here("code_repository/data/mc_burst_dyn.parquet")) ## dynamics data for mc burst
si_dyn_30.df <- read_parquet(here("code_repository/data/si_dyn_30.parquet")) ## single cue dynamics data (30 days)
dual_cue_dyn_30.df <- read_parquet(here("code_repository/data/dual_cue_dyn_30.parquet")) ## dual cue dynamics (30 days)

## validation data (comparing fitness to random spline)
validation.ls <- list.files(here("code_repository/data/si_validation/"), full.names = T)
validation.df <- do.call(rbind, lapply(validation.ls, read.csv))

# import in code
source(here("code_repository/functions/chabaudi_si_clean.R"))
source(here("code_repository/functions/par_to_df.R"))
source(here("code_repository/functions/chabaudi_si_clean_high.R"))
source(here("code_repository/functions/par_to_hm.R"))
source(here("code_repository/functions/par_to_hm_te.R"))

# color codes
orange <- "#fc8d59"
blue <- "#4575b4"

#=================================# # Cue perception mediates fitness #=================================# #———— A. conversion rate dynamics and fitness values ————# ## function for dynamics generation

# function for getting single cue infection dynamics
get_si_dyn <- function(df){
  ## processing model input
  par <- c(df$var1, df$var2, df$var3, df$var4) ## parameters
  cue <- df$cue ## cue choice
  log <- ifelse(df$log=="log", "log10", "none") ## log or not
  cue_range <- seq(df$low, df$high, by = df$by) ## cue range
  id <- df$id ## id
  
  ## get dynamics data
  dyn <- chabaudi_si_clean(
            parameters_cr = par, 
            parameters = parameters_tsukushi, 
            time_range = seq(0, 20, 0.001), 
            cue = cue, 
            cue_range = cue_range, 
            log_cue = log,
            immunity = "tsukushi",
            solver = "vode",
            dyn = TRUE)
  
  # append id
  dyn2 <- cbind(dyn, id = rep(id, nrow(dyn)))
  
  # return results
 return(dyn2)
}

run function to generate data

get 30 day simulation data (not used for fitness but for plotting purposes)

## get dynamics
si_opt.ls2 <- split(si_opt.df %>% filter(cue != "t"), seq(nrow(si_opt.df%>% filter(cue != "t"))))
si_dyn_30 <- mclapply(si_opt.ls2, get_si_dyn, mc.cores = 8)

## combine the dynamics file
si_dyn_30.df <- do.call(rbind, si_dyn_30)
si_dyn_30.df %>% group_by(id) %>% summarise(max_time = max(time))

## save
write_parquet(si_dyn_30.df, here("code_repository/data/si_dyn_30.parquet"))

process data to isolate conversion rate and pair them with fitness

si_opt.df[si_opt.df$cue == "t",]$fitness_20

[1] 9.7879

plot

## plot together
fitness_rank_cr.plt <- ggarrange(fitness_rank.plt, fitness_cr.plt, widths = c(1.2, 1), align = "h", common.legend = T,
                                 legend = "bottom")

Warning: Removed 11022 rows containing missing values or values outside the scale range (`geom_raster()`).
Warning: Removed 11022 rows containing missing values or values outside the scale range (`geom_raster()`).

#———— B. Reaction norms ————# ## function to obtain the cue values sensed by parasites from dynamics datas

get_rug <- function(df){
  ## process cue
  cue <- unique(df$cue)
  
  ## if cue contains "+", we need to first split them up and add them into the final dataframe
   if(stringr::str_detect(cue, "\\+")){
    cue_split <- stringr::str_split(string = cue, pattern = "\\+", simplify = T)
    ## get the two cues 
    cue_temp_1 <- cue_split[[1]]
    cue_temp_2 <- cue_split[[2]]
    ## filter dyn
    rug <- df %>% 
      filter(variable == cue_temp_1 | variable == cue_temp_2) %>% 
      dplyr::group_by(time) %>% 
      dplyr::mutate(sum = sum(value, na.rm = T)) %>% 
      select(time, value = sum, id)
  }
  
  # for cue with no addition, it is simply filtering for the values and returning it
  if(stringr::str_detect(cue, "\\+", negate = T)){
    rug <- df %>% 
      dplyr::filter(variable == cue) %>% 
      dplyr::select(time, value, id)}
  
  return(rug)
}

run function to obtain cue values

# join dynamics data with cue information
si_dyn_cue.df <- left_join(si_dyn.df, si_opt.df, by = "id")

# split based on individual labels
si_dyn_cue.ls <- si_dyn_cue.df %>% group_split(id)

# run function to get rug
si_rug <- mclapply(si_dyn_cue.ls, get_rug, mc.cores = 6)

# combine and save
si_rug.df <- do.call(rbind, si_rug)
write_parquet(si_rug.df, here("code_repository/data/si_rug.parquet"))

function to obtain reaction norm

get_si_rn <- function(df){
  ## read in the parameter sets
  par <- c(df$var1, df$var2, df$var3, df$var4)
  ## get cue range
  cue_range <- seq(df$low, df$high, by = df$by)
  ## nested function to convert parameter set into basis spline function
  rn <- par_to_df(par = par, cue_range = cue_range)
  
  # if parasite is sensing logged cue, we have to exponentiate it back so they are on the same scale! 
  rn2 <- rn
  if(stringr::str_detect(df$log, "log")){rn2$cue_range <- 10^(rn2$cue_range)}

  # append label to reaction norm, dyn, and rug
  rn2 <- data.frame(rn2, id = df$id)
  
  return(rn2)
}

run function to get reaction norms

# split dataframe of optimized rcues
si_opt.ls <- split(si_opt.df, seq(nrow(si_opt.df)))

# run function
si_rn <- lapply(si_opt.ls, get_si_rn)

# bind together
si_rn.df <- do.call(rbind, si_rn)
write_parquet(si_rn.df, here("code_repository/data/si_rn.parquet"))

preparing the dataset for plotting reaction norms and rugs

## we are only plotting relevent ranges of the reaction norm, meaning that these are the cue values actually "sensed" by the parasites in an infection. We can use the rug dataframe to get a rough estimate of that
si_rug_lim.df <- si_rug.df %>% 
  group_by(id)%>% 
  summarise(min_temp = min(value, na.rm = T)*0.9,
         max_temp = max(value, na.rm = T)*1.1) %>% ## calculate min and max for each cue
  left_join(si_opt.df, by = "id") %>% ## join with si_opt so we can get the cmmon cue
  ungroup() %>% 
  group_by(cue) %>% 
  summarise(min = min(min_temp), max = max(max_temp)) 

## process reaction norm data to restrict range and get label
si_rn.df_p <- si_rn.df %>% 
  left_join(select(si_opt.df, id, cue, log), by = "id") %>% ## get cue and log status from si_opt.df
  left_join(si_rug_lim.df, by = "cue") %>% ## via cue, get ranges we want to limit the rn to
  filter(cue_range <= max & cue_range >= min) %>% ## keep only cue values within range
  left_join(select(ez_label, id, cue_label), by = "id") # get labels


## refactor to reorder the order at which the cues are presented
si_rn.df_p$cue_label <- factor(si_rn.df_p$cue_label, 
                              levels = c("Asexual iRBC", "Sexual iRBC", "Total iRBC", "Gametocyte", "RBC"))

## split rug labels by none logged and logged
si_rug_cue.df <- si_rug.df %>% 
  distinct(id, value) %>% ## cut down on number of datapoints. For each id, keep only distinct points
  left_join(select(si_opt.df, id, cue, log), by = "id") %>% ## get cue and log status from si_opt.df
  left_join(select(ez_label, id, cue_label), by = "id") ## get cue label for plotting

## refactor rug
si_rug_cue.df$cue_label <- factor(si_rug_cue.df$cue_label, 
                              levels = c("Asexual iRBC", "Sexual iRBC", "Total iRBC", "Gametocyte", "RBC"))

### split
si_rug_cue.df_none <- si_rug_cue.df %>% filter(log == "none")
si_rug_cue.df_log <- si_rug_cue.df %>% filter(log == "log")

plot reaction norms and associated rug plots

## note that geom_point contains a set of scripts aimed at thinning the data points so they do not overcrowd!
fitness_rn.plt <- ggplot() +
  geom_line(data = si_rn.df_p, aes(x = cue_range, y = cr, color = log)) +
  geom_point(data = si_rn.df_p %>% 
  mutate(rel_cue = round(cue_range/(max-min)*100)) %>% 
    distinct(id, rel_cue, .keep_all = T) %>% 
    filter(rel_cue %% 10 ==1),
             aes(x = cue_range, y = cr, color = log, shape = log), size = 2) +
   geom_rug(data = si_rug_cue.df_none, aes(x = value), color = "#fc8d59", sides = "t", length = unit(0.1, "npc")) +
  geom_rug(data = si_rug_cue.df_log, aes(x = value), color = "#4575b4", sides = "b", length = unit(0.1, "npc")) +
  facet_wrap(~cue_label, scales = "free_x", ncol = 1) +
  labs(x = "Cue range", y = "Conversion rate", color = "Cue status", shape = "Cue status") +
  ylim(0, 1) +
  scale_color_manual(values=c( "#4575b4", "#fc8d59")) +
  scale_x_continuous(labels = function(x) format(x, scientific = T),
                     guide = guide_axis(check.overlap = TRUE)) + 
  theme_bw() +
  theme(axis.text.x = element_text(size = 8), legend.position = "bottom",
        
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank())

#———— Plotting fitness, dynamics, and rn together ————#

ggarrange(fitness_rank_cr.plt, fitness_rn.plt, widths = c(2.7, 1), align = "h", ncol = 2, labels = c("A", "B"))
ggsave(units = "px", dpi = 300, width = 2250, height = 1500, filename = here("code_repository/figures/fitness_rn.tiff"), bg = "white", scale = 1.25)

#=================================# # Dual cue fitness main figure #=================================# #————Relative ranking of dual cue vs single cue fitness———–# ## Prepare dual cue dataset

## For the dual cue models, we optimized each cue combinations in 2 ways, one with L-BFGS-B with 0.5X9 starting point and another with DEoptim + L-FGBS-B. For each cue combination, we are going to pick the strategy that gave us the higher fitness
dual_cue_f_final.df <- dual_cue_f_glb.df %>% 
  select(id, id_b, label, label_b, fitness_glb = fitness, 
         par1_glb = par1, par2_glb = par2, par3_glb = par3, par4_glb = par4, par5_glb = par5, par6_glb = par6,
         par7_glb = par7, par8_glb = par8, par9_glb = par9) %>% 
  left_join(
    select(dual_cue_f_lc.df,
         id, id_b, fitness_lc = fitness, 
         par1_lc = par1, par2_lc = par2, par3_lc = par3, par4_lc = par4, par5_lc = par5, par6_lc = par6,
         par7_lc = par7, par8_lc = par8, par9_lc = par9), by = c("id", "id_b")
  ) %>% ## rename columns and join global and local optimization fitness dataframes
  mutate(
    fitness = ifelse(fitness_glb > fitness_lc, fitness_glb, fitness_lc),
    par1 = ifelse(fitness_glb > fitness_lc, par1_glb, par1_lc),
    par2 = ifelse(fitness_glb > fitness_lc, par2_glb, par2_lc),
    par3 = ifelse(fitness_glb > fitness_lc, par3_glb, par3_lc),
    par4 = ifelse(fitness_glb > fitness_lc, par4_glb, par4_lc),
    par5 = ifelse(fitness_glb > fitness_lc, par5_glb, par5_lc),
    par6 = ifelse(fitness_glb > fitness_lc, par6_glb, par6_lc),
    par7 = ifelse(fitness_glb > fitness_lc, par7_glb, par7_lc),
    par8 = ifelse(fitness_glb > fitness_lc, par8_glb, par8_lc),
    par9 = ifelse(fitness_glb > fitness_lc, par9_glb, par9_lc)
  ) ## the columns are assigned based on which strategy produced the highest fitness

## write csv
write.csv(dual_cue_f_final.df, here("code_repository/data/dual_cue_fitness_final.csv"))

Combine single cue and dual cue dataset

dual_cue_f_final.df <- read.csv(here("code_repository/data/dual_cue_fitness_final.csv"))

dual_si_fitness.df <- dual_cue_f_final.df %>% 
  select(id, id_b, label, label_b, fitness_dual = fitness) %>% 
  left_join(select(si_opt.df, id, fitness_si = fitness_20), by = "id") %>% ## get fitness of single cue model based on first cue
  left_join(select(si_opt.df, id_b = id, fitness_si_b = fitness_20), by = "id_b") %>%  ## get fitness of single cue model based on second cue
  mutate(fitness_si_final = ifelse(fitness_si > fitness_si_b, fitness_si, fitness_si_b)/9.883602,
         label_comb = paste(label, "&", label_b), ## select the highest single cue model fitness
         fitness_dual_norm = fitness_dual/9.883602)## normalize by time fitness

plot fitness rank

#————Time series conversion rate of dual cue models———–# Note that our dual cue models take 9 parameters. To ensure that any difference between dual and single cue models is due to the inclusion of additional cues and NOT higher spline flexibility,we performed L-BFGS-B (local) optimization of single cue models with df = 9 to control for the spline flexibility.

dynamics simulation of high parameter cues

## best dual cue model: I log and R log
Rlog_Ilog.cr <- chabaudi_si_clean(
  parameters_cr = c(4.446192033,    10.97518275,    1.38762817, 23.3059254, -3.452052371,   -18.0070692,    39.66614226,    -3.545193141,   18.78350799),
  immunity = "tsukushi",
  parameters = parameters_tsukushi,
  time_range = seq(0, 20, by = 1e-3),
  cue_range =  seq(6, 7, by = 1/500),
  cue_range_b = seq(0, log10(6*(10^6)), by = (log10(6*(10^6)))/500),
  cue = "R",
  cue_b = "I",
  log_cue = "log10",
  log_cue_b = "log10",
  solver = "vode",
  dyn = T
)

## when time is used as a cue (9 parameters). chabaudi_si_clean_high is simply a variation of chabaudi_si_clean that makes a more flexible spline!
time_high.cr <- chabaudi_si_clean_high(
  parameters_cr = c(9.154314,  -7.570829, -22.506638 ,  3.382405 ,-13.453519 ,-17.011485  , 3.678181, -12.851895 ,-26.115158),
  immunity = "tsukushi",
  parameters = parameters_tsukushi,
  time_range = seq(0, 20, by = 1e-3),
  cue_range =  seq(0, 20, by = 1e-3),
  cue = "t",
  solver = "vode",
  dyn = T)

## when asexual iRBC is used as a cue (high flexibility)
I_high.cr <- chabaudi_si_clean_high(
  parameters_cr = c(1.296675,  3.544034 , 4.907484,  2.174249, -3.238309 ,-5.181614 ,-1.645072 , 1.834302 , 1.581011),
  immunity = "tsukushi",
  parameters = parameters_tsukushi,
  time_range = seq(0, 20, by = 1e-3),
  cue_range =  seq(0, log10(6*(10^6)), by = (log10(6*(10^6)))/5000),
  cue = "I",
  log_cue = "log10",
  solver = "vode",
  dyn = T)

## when RBC is used as cue (high flexibility)
R_high.cr <- chabaudi_si_clean_high(
  parameters_cr = c(5.0340348 ,  0.5846168 ,  0.3749648 ,  0.6842673  , 2.4748107 , 10.9036034 , 16.8246316, -24.8690971 ,1.8007238),
  immunity = "tsukushi",
  parameters = parameters_tsukushi,
  time_range = seq(0, 20, by = 1e-3),
  cue_range =  seq(log10(10^6), log10(10^7), by = (log10(10^7)-log10(10^6))/5000),
  cue = "R",
  log_cue = "log10",
  solver = "vode",
  dyn = T)

plot

dual_selected_cr.pl <- ggplot() +
  geom_line(data = dual_selected_cr.df, aes(color = label_new, x = time, y = value), size = 1) +
  geom_point(data = dual_selected_cr.df %>% filter(time%%1 == 0), aes(color = label_new, x = time, y = value, shape = label_new), size = 3) +
  labs(x = "Time (days)", y = "Conversion rate", color = "Cue(s)", shape = "Cue(s)") +
  xlim(0, 20) +
  scale_color_manual(values = c("#fc8d59","#fdcb44","black", "#4575b4")) +
  theme_classic() +
  theme(legend.position="right",
        plot.margin = margin(t = 40, r = 0, b = 0, l = 0, unit = "pt")) +
  guides(color = guide_legend(nrow = 4, byrow = TRUE))

Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.
This warning is displayed once every 8 hours.
Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.

#——– Reaction norm heatmap of R log10 + I log10 ————# # Process data

# make heatmap data.frame
Rlog_Ilog.hm <- par_to_hm_te(par = c(4.446192033,   10.97518275,    1.38762817, 23.3059254, -3.452052371,   -18.0070692,    39.66614226,    -3.545193141,   18.78350799),
             cue_range = seq(6, 7, length.out = 500),
             cue_range_b = seq(0,   6.77815125, length.out = 500))

# process dynamics
Rlog_Ilog.dyn <- Rlog_Ilog.cr %>% 
  tidyr::pivot_wider(names_from = variable, values_from = value) %>% 
  mutate(log_R = log10(R),
         log_I = log10(I))

plot

Rlog_Ilog_rn.pl <- ggplot() +
  geom_raster(data = Rlog_Ilog.hm, aes(x = cue_range_b, y = cue_range, fill = cr)) +
  scale_fill_viridis_c() +
  geom_path(data = Rlog_Ilog.dyn, aes(x = log_I, y = log_R), color = "white", arrow = arrow(angle = 30, length = unit(0.1, "inches"))) +
  geom_point(data = Rlog_Ilog.dyn %>% filter(row_number() %% 1000 == 1 & time <= 20), aes(x = log_I, y = log_R), color = "white") +
  xlim(0.99*min(hablar::s(Rlog_Ilog.dyn$log_I), na.rm = T), 1.01* max(hablar::s(Rlog_Ilog.dyn$log_I), na.rm = T)) +
  ylim(0.99*min(hablar::s(Rlog_Ilog.dyn$log_R), na.rm = T),1.01* max(hablar::s(Rlog_Ilog.dyn$log_R), na.rm = T)) +
  labs(y = "RBC log10", x = "Asexual iRBC log10", fill = "Conversion\nrate") +
  theme_classic() +
  theme(legend.position = "right")

#———– Plot together ————#

#=================================# # Dual cue conversion rate supplementary figure #=================================# ## Function to simulate dual cue dynamics

dual_cue_dyn <- function(df){
  ## process cues
  cue <- df$cue
  cue_b <- df$cue_b
  
  ## process log
  log <- ifelse(str_detect(df$id, "log"), "log10", "none")
  log_b <- ifelse(str_detect(df$id_b, "log"), "log10", "none")
  
  # process cue_range. ensure that both cue ranges are of the same length
  cue_range <- seq(df$low, df$high, length.out = 500)
  cue_range_b <- seq(df$low_b, df$high_b, length.out = 500)
  
  # get parameter set
  par <- c(df$par1, df$par2, df$par3, df$par4, df$par5, df$par6, df$par7, df$par8, df$par9)
  
  # simulate dynamics
  dyn <- chabaudi_si_clean(
    parameters_cr = par,
    immunity = "tsukushi",
    parameters = parameters_tsukushi,
    time_range = seq(0, 20, 0.01),
    cue = cue,
    cue_b = cue_b,
    cue_range = cue_range,
    cue_range_b = cue_range_b,
    log_cue = log,
    log_cue_b = log_b,
    solver = "vode",
    gam = "te",
    dyn = T)
  
  # return 
  dyn2 <- cbind(id = df$id, id_b = df$id_b, 
                label = df$label, label_b = df$label_b,
                cue = cue, cue_b = cue_b, dyn)
  
  write_parquet(dyn2, here(paste0("code_repository/data/dual_cue_dyn/", df$id, "_", df$id_b, ".parquet")))
}

Get dynamics of dual cue models

## For the dual cue fitness dataframe, add the cue range for each individual cue
dual_cue_f_final.df_p <- dual_cue_f_final.df %>% 
  left_join(select(cue_range_si_alt.df, id, cue, low, high), by = "id") %>% 
  left_join(select(cue_range_si_alt.df, id_b = id, cue_b = cue, low_b = low, high_b = high), by = "id_b")

## Split dataframes
dual_cue_f_final.ls <- split(dual_cue_f_final.df_p, seq(nrow(dual_cue_f_final.df_p)))

## Run function
mclapply(dual_cue_f_final.ls, dual_cue_dyn, mc.cores = 6)

## Concat all files
dual_cue_dyn.ls <- list.files(path = here("code_repository/data/dual_cue_dyn"), pattern = "*.parquet", full.names = T)
dual_cue_dyn.df <- do.call(rbind, lapply(dual_cue_dyn.ls, read_parquet))
write_parquet(dual_cue_dyn.df, here("code_repository/data/dual_cue_dyn.parquet"))

get 30 days dynamic (for plotting purposes)

## Run function
mclapply(dual_cue_f_final.ls, dual_cue_dyn, mc.cores = 6)

## Concat all files
dual_cue_dyn_30.ls <- list.files(path = here("code_repository/data/dual_cue_dyn_30"), pattern = "*.parquet", full.names = T)
dual_cue_dyn_30.df <- do.call(rbind, lapply(dual_cue_dyn_30.ls, read_parquet))
write_parquet(dual_cue_dyn_30.df, here("code_repository/data/dual_cue_dyn_30.parquet"))

Preparing dataset for plotting

## filter out only conversion rate and attach fitness values
dual_cue_cr.df <- dual_cue_dyn.df %>% 
  filter(variable == "cr") %>% 
  left_join(select(dual_cue_f_final.df, id, id_b, fitness), by = c("id", "id_b")) %>% 
  mutate(label_comb = paste(label, "&", label_b))

## Sanity check that the dynamics produced the same fitness as the optimized values. Yes!
dual_cue_dyn.df %>% 
  filter(variable == "tau_cum") %>% 
  filter(time == 20) %>% 
  distinct(id, id_b, value) %>% 
  left_join(select(dual_cue_f_final.df, id, id_b, fitness), by = c("id", "id_b")) %>% 
  mutate(diff = value-fitness)

Plot

ggsave(here("code_repository/figures/dual_cue_cr.tiff"), units = "px", width = 2000, height = 1500, dpi=300,  bg = "white")

Warning: Removed 4080 rows containing missing values (`geom_raster()`).

#=================================# # Experimental disease maps of P. chabaudi #=================================# ## Import in data

# import in https://academic.oup.com/emph/article/2018/1/127/5045871?login=true
## (2018 published in EMPH)
emph_2018 <- readxl::read_xls(here("code_repository/experimental_data/Huijben_2018_EMPH.xls"), sheet = 1)

# import in https://onlinelibrary.wiley.com/doi/10.1111/j.1558-5646.2010.01068.x
## (2010 published in Evolution)
evo_2010 <- readxl::read_xls(here("code_repository/experimental_data/Huijben_2010_evolution.xls"), sheet = 1)

# import in https://journals.plos.org/plospathogens/article?id=10.1371/journal.ppat.1003578#:~:text=The%20philosophy%20is%20that%20aggressive,longer%20feel%20sick%20%5B13%5D.
## (2013 in PLoS pathogen)
plos_2013_1 <- readxl::read_xlsx(here("code_repository/experimental_data/Huijben_2013_PLoS.xlsx"), sheet = 2)

plos_2013_2 <- readxl::read_xlsx(here("code_repository/experimental_data/Huijben_2013_PLoS.xlsx"), sheet = 3)

# import in https://onlinelibrary.wiley.com/doi/10.1111/j.1420-9101.2011.02369.x
## (2011 in Journal of Evolutionary Biology)
eseb_2011 <- readxl::read_xls(here("code_repository/experimental_data/Huijben_2011_eseb.xls"), sheet = 1)

# import in https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3939351/
## (2011 in Journal of American naturalist). Note private dataset so not provided in the supplementary!
amna_2011 <- readxl::read_xls(here("experimental_data/Pollitt_2011_naturalist.xls"), sheet = 1)

Clean data

## for EMPH 2018 study, include only infection series without drugs were R-inoculum is administered by itself, which includes 6, 7, 8, 9, 10. Box 6 has a starting inoculum number of 10^6, which is most similar to other studies. Filtering between day 3-21 because those are the days where we have single day data.
emph_2018_ss.df <- emph_2018 %>% 
  filter(Box %in% seq(6, 10) &
         dplyr::between(Day, 3, 21)) %>% 
  mutate(dose = case_when(
    Box == 6 ~ 10^6,
    Box == 7 ~ 10^5,
    Box == 8 ~ 10^3,
    Box == 9 | Box == 10 ~ 10^1
  )) %>% 
  mutate(strain = "As6p",
         study = "emph2018",
         study_strain = paste0(strain, study),
         id = paste0(study, strain, Box, Mouse, 1),
         RBC = RBC * (10^6)) %>%
  select(day = Day,
         mouse = Mouse, 
         RBC, 
         asex = Rasex,
         gam = Rgam,
         dose,
         strain,
         study,
         study_strain,
         id)

## for 2011 eseb, only day 3-17 data are analyzed because those are the days where gametocyte data are available
eseb_2011_ss.df <- eseb_2011 %>% 
  filter(Clones == "R" & between(Day, 3, 17) &
           Drugs == "N") %>% 
  mutate(dose = 10^6,
         strain = "As8p",
         study = "eseb2011",
         study_strain = paste0(strain, study),
         RBC = RBC*(10^6),
         id = paste0(study, strain, Box, Mouse, 2)) %>% 
  select(day = Day,
         mouse = Mouse,
         RBC,
         asex = R.asex,
         gam = R.gam,
         dose,
         strain,
         study,
         study_strain,
         id)

## for evolution_2010, single infection data for both resistant and susceptible clones are available without drug treatment
evolution_2010_ss.df <- evo_2010 %>% 
  filter(Clone == "R" | Clone == "S") %>%
  filter(between(Day, 3, 21) &
           Drugs == "nodrugs") %>% 
  mutate(asex = R.asex + S.asex,
         gam = R.gam + S.gam,
         dose = 10^6,
         RBC = RBC*(10^6),
         study = "evol2011",
         strain = ifelse(Clone == "R", "As12", "AJ51"),
         study_strain = paste0(strain, "_", study),
         id = paste0(study, strain, Box, Mouse, 3)) %>% 
  select(day = Day,
         mouse = Mouse,
         RBC,
         asex,
         gam,
         dose,
         strain,
         study,
         study_strain,
         id)

## for amnat 2011, get single infection data. Filter out any mice that have missing data. Set negative asexuasl data ot 0
amna_2011_ss.df <-  amna_2011 %>% 
  filter(treat %in% c("AJ", "AS", "ER", "CR", "CW", "DK")) %>% 
  mutate(asex = tot.para - tot.gcyte,
         gam = tot.gcyte,
         dose = 10^6,
         study = "amna_2011",
         RBC = rbc/(10^6),
         study_strain = paste0(treat, "_", study),
         id = paste0(study, treat, div, mouse, 4)) %>% 
  mutate(asex = ifelse(asex < 0, 0, asex)) # sometimes total parasite is less than gametocyte so need to correct for this

### check for NA by groups
amna_na.id <- amna_2011_ss.df %>% 
  filter_at(vars(asex, gam, RBC), all_vars(is.na(.))) %>% 
  distinct(id) %>% 
  select(id)

amna_2011_ss.df2 <- amna_2011_ss.df %>% 
  filter(!(id %in% amna_na.id$id)) %>% 
  select(day,
         mouse,
         RBC,
         asex,
         gam,
         dose,
         strain = treat,
         study,
         study_strain,
         id)

## rbind
exp_ss.df <- rbind(emph_2018_ss.df, eseb_2011_ss.df, evolution_2010_ss.df, amna_2011_ss.df2)

## write
write.csv(exp_ss.df, here("code_repository/data/experimental_data.csv"))

Prepare dataset for plotting

names(exp_ss.df) <- c("X", "Day", "Mouse", "RBC", "iRBC", "Gametocyte", "Dose", "Strain", "Study", "Study_strain", "id")

## prepare a list of variable combinations we want to plot
exp_var.comb <- tidyr::expand_grid(x = c("RBC", "iRBC", "Gametocyte"),
                   y = c("RBC", "iRBC", "Gametocyte")) %>% ## get all pairwise combinations of variables
  filter(x != y) %>% ## remove incidences where the 2 variables are the same
  mutate(tmp = paste0(pmin(x, y), pmax(x, y))) %>% ## eliminate same variable but different order
  slice_head(n = 1, by = tmp) %>% 
  select(-tmp) 

List-wise plotting

## x and y axis are not logged!
exp_xy.pl_ls <- map2(exp_var.comb$x, exp_var.comb$y, ~ {
  x_col <- .x
  y_col <- .y
  
  ggplot(exp_ss.df, aes_string(x = x_col, y = y_col)) +
        geom_path(aes(colour = Day, group = id), arrow = arrow(type = "closed", angle = 10, length = unit(0, "inches"))) +
        theme_classic() + 
        scale_color_viridis_c(option = "A", limits = c(3, 21)) +
        labs(color = "Days post-infection")  +
        scale_y_continuous(labels = function(x) format(x, scientific = TRUE)) +
        scale_x_continuous(labels = label_scientific(digits = 1))
}
)

## x-axis is logged
exp_xlogy.pl_ls <- map2(exp_var.comb$x, exp_var.comb$y, ~ {
  x_col <- .x
  y_col <- .y
  
  ggplot(exp_ss.df, aes_string(x = sprintf("log10(%s)", x_col), y = y_col)) +
        geom_path(aes(colour = Day, group = id), arrow = arrow(type = "closed", angle = 10, length = unit(0, "inches"))) +
        theme_classic() + 
        scale_color_viridis_c(option = "A", limits = c(3, 21)) +
        labs(color = "Days post-infection", x = paste(x_col, "log"))  +
        scale_y_continuous(labels = function(x) format(x, scientific = TRUE))
}
)

## y-axis logged
exp_xylog.pl_ls <- map2(exp_var.comb$x, exp_var.comb$y, ~ {
  x_col <- .x
  y_col <- .y
  
  ggplot(exp_ss.df, aes_string(x = x_col, y = sprintf("log10(%s)", y_col))) +
        geom_path(aes(colour = Day, group = id), arrow = arrow(type = "closed", angle = 10, length = unit(0, "inches"))) +
        theme_classic() + 
        scale_color_viridis_c(option = "A", limits = c(3, 21)) +
        labs(color = "Days post-infection", y = paste(y_col, "log"))   +
        scale_x_continuous(labels = label_scientific(digits = 1))
}
)

## both x and y axis is logged
exp_xlogylog.pl_ls <- map2(exp_var.comb$x, exp_var.comb$y, ~ {
  x_col <- .x
  y_col <- .y
  
  ggplot(exp_ss.df, aes_string(x = sprintf("log10(%s)", x_col), y = sprintf("log10(%s)", y_col))) +
        geom_path(aes(colour = Day, group = id), arrow = arrow(type = "closed", angle = 10, length = unit(0, "inches"))) +
        theme_classic() + 
        scale_color_viridis_c(option = "A", limits = c(3, 21)) +
        labs(color = "Days post-infection", x = paste(x_col, "log"), y = paste(y_col, "log")) 
}
)

## plot together
ggarrange(plotlist = c(exp_xy.pl_ls, exp_xlogy.pl_ls, exp_xylog.pl_ls, exp_xlogylog.pl_ls), 
          common.legend = T, align = "hv")
ggsave(here("code_repository/figures/exp_disease-curve.tiff"), units = "px", width = 2250, height = 1500, scale = 1.4, dpi=300,  bg = "white")

#=================================# # Simulated disease curve graph #=================================# ## Function to obtain dynamics data based on the dual cue input

get_dual_rn <- function(df){
  ## assign the two cues
  cue <- unique(df$cue)
  cue_b <- unique(df$cue_b)
  
  ## assign log status
  log <- ifelse(str_detect(unique(df$id), "log"), "log", "none")
  log_b <- ifelse(str_detect(unique(df$id_b), "log"), "log", "none")
  
  ## assign which cues are going to be logged
  if(log == "log" & log_b == "none"){logged_cue <- cue}
  if(log == "none" & log_b == "log"){logged_cue <- cue_b}
  if(log == "log" & log_b == "log"){logged_cue <- c(cue, cue_b)}
  if(log == "none" & log_b == "none"){logged_cue <- c()}

  ## keep variables that corresponds to the cue used
  ### for dataframes that does not involve combined variables such as I+Ig
  if(isTRUE(str_detect(cue, "\\+", negate = T)) & isTRUE(str_detect(cue_b, "\\+", negate = T))){
    df_f <- df %>% 
      filter(variable %in% c(cue, cue_b, "cr")) %>% 
      mutate(value = case_when(
        variable %in% logged_cue ~ log10(value),
        TRUE ~ value
      )) %>% ## log transform values only when they match with the logged cue list
      filter(value >= 0) ## filter out values <0, these happen due to stiffness of models but are not relevant
         
  } else{
    ### assign both cues to a list
    cue_ls <- c(cue, cue_b)
    ### pick the cue that has the "+" sign
    combined_cue <- cue_ls[grepl("\\+", cue_ls)]
    non_combined_cue <- cue_ls[!grepl("\\+", cue_ls)] ### this is the none combined cue
    ### unlist the cues
    cue_unlist <- unlist(str_split(combined_cue, "\\+"))
    
    ### get filtered dataset containing only non-combine cue
    df_f1 <- df %>% 
      filter(variable %in% c(non_combined_cue, "cr"))
    
    ### get filtered dataset containing combined cue. These will be summed up and bound back to the previous
    df_f2 <- df %>% 
      filter(variable %in% cue_unlist) %>% ## keep only variables that we will combine
      group_by(time) %>% ## for each time point, group the variables
      mutate(value = sum(value, na.rm = T),
             variable = combined_cue) %>% ## recalculate the value as sum of the values and reassign variable!
      distinct(time, .keep_all = T) ## note that because we are mutating we must dedeuplicate the records
    
    #### combine the two and log transform if necessary
    df_f <- rbind(df_f1, df_f2) %>% 
      mutate(value = ifelse(variable %in% logged_cue, log10(value), value)) %>% ## log transform values only when they match with the logged cue list
      filter(value >= 0) 
  }
  
  ## Convert dataframes wider such that the different variables have their own columns
  df_fp <- df_f %>% 
    mutate(variable_id = ifelse(variable == cue, paste0(variable, "_", log), paste0(variable, "_", log_b))) %>%  ### assign a unique variable id that could later be used to assign labels
    left_join(select(ez_label, id, long_label), by = c("variable_id" = "id")) %>% 
    mutate(long_label = ifelse(variable == "cr", "cr", long_label)) %>% ## manually add cr
    mutate(long_label = gsub(" ", "_", long_label)) %>%  ## convert spaces to _ for plotting
    pivot_wider(names_from = long_label, values_from = value, id_cols = c(time, id, id_b)) %>% 
    filter(time >= 1) %>% ## filter out day 0->1 because all cr = 0 before that
    arrange(time) ## this is needed to prevent geom_path from joining the first and last data point
  
  ## assign NAs (meaning no stuff is produced yet to 0)
  df_fp[is.na(df_fp)] <- 0

  return(df_fp)
}

Run function to get a curated dataset containing only relevant

## split dual dynamics dataframe into list grouped by the dual cues
dual_cue_dyn.ls <- dual_cue_dyn.df %>% group_split(id, id_b)

## run function across list
dual_cue_rn.ls <- mclapply(dual_cue_dyn.ls, get_dual_rn, mc.cores = 6)

## sanity checks that we are actually summing iRBCs. Rings out
max(dual_cue_rn.ls[[2]]$Total_iRBC)
max(dual_cue_rn.ls[[4]]$Asexual_iRBC)
max(dual_cue_rn.ls[[6]]$Sexual_iRBC)

max(dual_cue_rn.ls[[10]]$Total_iRBC)
max(dual_cue_rn.ls[[12]]$Asexual_iRBC)
max(dual_cue_rn.ls[[14]]$Sexual_iRBC)

plot

## list apply all dataframes
cue_cue_rn_pl.ls <- lapply(dual_cue_rn.ls,
       function(x){
          ## get names of columns used in the x and y axis
         axis_cols <- setdiff(names(x), c("time", "id", "id_b", "cr"))
         
         ## ggplot
         ggplot() +
           geom_path(data = x, aes_string(x = axis_cols[[1]], y = axis_cols[[2]], color = "cr"),
                     arrow = arrow(length = unit(c(rep(0, nrow(x) - 2), 0.25), "inches")),
                     size = 1.5) +
           geom_point(data = x %>% filter(time %% 1 == 0), 
                      aes_string(x = axis_cols[[1]], y = axis_cols[[2]]), size = 1.5, shape = 1) +
           theme_classic() +
           scale_color_viridis_c(limits = c(0, 1)) +
           labs(color = "Conversion rate", x = gsub("_", " ", axis_cols[[1]]), y = gsub("_", " ", axis_cols[[2]])) +
           scale_x_continuous(labels = label_scientific(digits = 2)) +
           scale_y_continuous(labels = label_scientific(digits = 2))
       })

## arrange together
cue_cue_rn.pl <- ggarrange(plotlist = cue_cue_rn_pl.ls, ncol = 5, nrow = 8, common.legend = T, align = "hv") 
ggsave(here("code_repository/figures/sim_disease-curve.tiff"), units = "px", width = 2250, height = 2500, scale = 2.2, dpi=300,  bg = "white")

#=================================# # Curating list of selected exo and simulated # disease curves (main figure) #=================================# # plot

# main figure
ggarrange(
   cue_cue_rn_pl.ls[[31]] +
    scale_x_continuous(labels = scales::number_format(accuracy = 0.1)) +
    scale_y_continuous(labels = scales::number_format(accuracy = 1)) +
    theme(legend.position = "none"), 
          
   exp_xlogylog.pl_ls[[1]] + theme(legend.position = "none"),

   cue_cue_rn_pl.ls[[35]] + 
     scale_y_continuous(labels = scales::number_format(accuracy = 1)) +
     theme(legend.position = "none",
           axis.text.x=element_text(size=6.5)), 
   
   exp_xylog.pl_ls[[1]] + 
     scale_y_continuous(labels = scales::number_format(accuracy = 1)) +
     theme(legend.position = "none",
           axis.text.x=element_text(size=6.5)),
   
   cue_cue_rn_pl.ls[[30]] + 
     scale_x_continuous(labels = scales::number_format(accuracy = 0.1)) +
     theme(legend.position = "none"), 
   
   exp_xlogy.pl_ls[[1]] + theme(legend.position = "none"), 
   
   cue_cue_rn_pl.ls[[15]] + 
     scale_x_continuous(labels = scales::number_format(accuracy = 0.1)) +
     theme(legend.position = "none"), 
   
   exp_xlogy.pl_ls[[2]] + theme(legend.position = "none"),
   
   align = "hv", ncol = 4, nrow = 2
)

Scale for x is already present.
Adding another scale for x, which will replace the existing scale.
Scale for y is already present.
Adding another scale for y, which will replace the existing scale.
Scale for y is already present.
Adding another scale for y, which will replace the existing scale.
Scale for x is already present.
Adding another scale for x, which will replace the existing scale.
Scale for x is already present.
Adding another scale for x, which will replace the existing scale.

ggsave(here("code_repository/figures/sim_exp_disease_curve_main.tiff"), units = "px", width = 2250, height = 1000, scale = 1.25, dpi=300,  bg = "white")

# get legend separately
ggarrange(
   cue_cue_rn_pl.ls[[31]] +
    scale_x_continuous(labels = scales::number_format(accuracy = 0.1)) +
    scale_y_continuous(labels = scales::number_format(accuracy = 1)) +
    theme(legend.position = "top"), 
          
   exp_xlogylog.pl_ls[[1]] + theme(legend.position = "top"))

Scale for x is already present.
Adding another scale for x, which will replace the existing scale.
Scale for y is already present.
Adding another scale for y, which will replace the existing scale.

ggsave(here("code_repository/figures/sim_exp_disease_curve_legend.tiff"), units = "px", width = 2250, height = 500, scale = 1.25, dpi=300,  bg = "white")

#=================================# # MC simulation of single cue and dual cue infection #=================================# #———— Impact of all parameter variation on fitness ————# ## Append all fitness data and sanity checks

## Import all fitness data and make into single dataframe
mc_all.ls <- list.files(path = here("code_repository/data/mc_all_fitness"), pattern = "*.csv", full.names = T)

## filter out MC records containing R log & I log. These data have different headers and would need to be processed differently
mc_all_sc.ls <- mc_all.ls[!grepl("R_log10-I_log10*", mc_all.ls)]
mc_all_Rlog_Ilog.ls <- mc_all.ls[grepl("R_log10-I_log10*", mc_all.ls)]

## we are expecting 11*5010 = 55110 data files
length(mc_all_sc.ls)
length(mc_all_Rlog_Ilog.ls) ## 5010 records nice

## read and append
mc_all_fitness_rc.df <- do.call(rbind, mclapply(mc_all_sc.ls, function(x){df <- read.csv(x)}, mc.cores = 6))
mc_all_fitness_Rlog_Ilog.df <- do.call(rbind, mclapply(mc_all_Rlog_Ilog.ls, function(x){df <- read.csv(x)}, mc.cores = 6))

## just correct the log -> log10 for consistency purpose and also rename id to iter
mc_all_fitness.df <- mc_all_fitness_Rlog_Ilog.df %>% 
  mutate(cue = "R log & I log", log = "log10") %>% ## renaming the cue names of R log and I log data so they fit
  select(-c(cue_b, log_b)) %>% 
  rbind(mc_all_fitness_rc.df) %>% 
  mutate(log = ifelse(log == "log10", "log", "none"),
         iter = id) %>% select(-id)

## write
write_parquet(mc_all_fitness.df, here("code_repository/data/mc_all_fitness.parquet"))

## for each cue and log, id = 1 is where all parameters are default values. Given that we simulated the time step with 0.01 rather than  0.001, we need to check whether this higher time step alters the fitness values. We can see that the difference is very small ~0.002 so this will not affect any our outcomes.
mc_all_fitness.df %>% filter(iter == 1) %>% 
  left_join(select(si_opt.df, cue, log, fitness_20), by = c("cue", "log")) %>% 
  mutate(diff = max_fitness - fitness_20)

process data for plotting fitness

## attach label
mc_all_fitness.df_p <- mc_all_fitness.df %>% 
  left_join(ez_label, by = c("cue", "log")) %>% # get label
  mutate(long_label = ifelse(cue == "R log & I log", "RBC log &\nasexual iRBC log", long_label),
         long_label = ifelse(cue == "t", "Time", long_label)) ## assign label manually to dual cue data

## manually check that everything is assigned correctly
mc_all_fitness.df_p %>% distinct(long_label, cue, log)

## get the reference fitness (default parameter variation), which is where iter = 1
mc_all_fitness_ref.df <- mc_all_fitness.df_p %>% filter(iter == 1)

## get the rest of the data points (excluding iter == 1) and calculate median and mean
mc_all_fitness_rand.df <- mc_all_fitness.df_p %>% 
  filter(iter != 1) 

## get mean and mode in a separate df
mc_all_fitness_sum.df <- mc_all_fitness_rand.df %>% 
  group_by(long_label) %>% 
  summarize(mean = mean(max_fitness),
            median = median(max_fitness),
            geom_mean = exp(mean(log(max_fitness)))) 

plot fitness variation

#———— Impact of individual parameter variation on fitness ————# ## combine all single parameter variation files

## list of file paths linked to the single parameter files
mc_single.ls <- list.files(path = here("code_repository/data/mc_single_fitness"), pattern = "*.csv", full.names = T)

## filter out R log and I log data. these will be attached later
mc_single_sc.ls <- mc_single.ls[!grepl("R_log10-I_log10*", mc_single.ls)]
mc_single_Rlog_Ilog.ls <- mc_single.ls[grepl("R_log10-I_log10*", mc_single.ls)]

## check number of files
length(mc_single_sc.ls)
length(mc_single_Rlog_Ilog.ls)

## read and combine
mc_single_sc_fitness.df <- do.call(rbind, mclapply(mc_single_sc.ls, function(x) read.csv(x), mc.cores = 6))
mc_single_Rlog_Ilog_fitness.df <- do.call(rbind, mclapply(mc_single_Rlog_Ilog.ls, function(x) read.csv(x), mc.cores = 6))

## change id -> iter and correcting log label
mc_single_fitness.df <- mc_single_Rlog_Ilog_fitness.df %>% 
  mutate(cue = "R log & I log") %>% ## manually recode cue for the dual cue results
  select(-c(cue_b, log_b)) %>% ## these columns are not present in the single cue models and are removed
  rbind(mc_single_sc_fitness.df) %>% 
  mutate(log = ifelse(log == "log10", "log", "none"),
                                                  iter = id) %>% select(-id)
## write
write_parquet(mc_single_fitness.df, here("code_repository/data/mc_single_fitness.parquet"))

## sanity check. iter = 1 is where all parameters are default. Should produce the same fitness values!
mc_single_fitness.df %>% filter(iter == 1)

## sanity check. We should have the same values for iter 1 across all and single parameter datasets. 
mc_single_fitness.df %>% filter(iter == 1)

process data for plotting

## note that for each iter across the "single" and "all" dataset, the parameter alteration is the same. Thus, for each data point at which all parameter are varied, there is a corresponding datapoint where only one parameter is varied. We can join these 2 dataset by iter and cue and log
mc_single_all_fitness.df <- mc_single_fitness.df %>% 
  left_join(select(mc_all_fitness.df, fitness_all = max_fitness, cue, log, iter), by = c("cue", "log", "iter"))

## make the dataframe into a long format such that all fitness variations are in a single column
mc_single_fitness.long <- mc_single_all_fitness.df %>% 
  filter(iter != 1) %>% ## filter out iter = 1, which does not have variation
  left_join(select(mc_all_fitness_ref.df, fitness_ref = max_fitness, cue , log),
            by = c("cue", "log")) %>% ## add in the deterministic fitness values
  select(-c("X", "rho", "burst", "iota_N1", "iota_N2", "phi_N1", "phi_N2")) %>% # keep only fitness and associated labels
  tidyr::pivot_longer(-c("cue", "log", "fitness_all", "iter", "fitness_ref")) %>% ## make long
  mutate(parameter = gsub("fitness_", "", name)) ## isolate parameter being altered

## calculate degree deviation from deterministic values. note that rel_diff_single ranges from -1 to 1. -1 ->  one variable perturbation is acting in the opposite direction to the overall perturbation caused by randomizing all variables. 0 parameter variation contributes very little, 1 -> one variable contributes a lot
mc_single_fitness.long_p <- mc_single_fitness.long %>% 
  mutate(diff_single = value-fitness_ref, ## pertubation to fitness caused by single parameter variation
         diff_all = fitness_all-fitness_ref, ## pertubation to fintess caused by all parameter varying
         rel_diff_single = diff_single/diff_all ## normalized pertunation to fitness (single parameter)
         ) %>% 
  left_join(ez_label, by = c("cue", "log")) %>% # get label
  mutate(long_label = ifelse(cue == "R log & I log", "RBC log &\nasexual iRBC log", long_label)) ## maually assign label to dual cue data

## calculate summary statistics. this includes median, credible interval (contains 89% data points calculated via Highest Density Interval, which is better for skewed data)
mc_single_fitness.sum <- mc_single_fitness.long_p %>% 
  group_by(long_label, parameter) %>% 
  summarise(ci_lower = ci(rel_diff_single, method = "HDI", ci = 0.89)[[2]],
            ci_higher = ci(rel_diff_single, method = "HDI", ci = 0.89)[[3]],
            quantile_low = quantile(rel_diff_single, 0.025),
            quantile_high = quantile(rel_diff_single, 0.975),
            median = median(rel_diff_single),
            mean = mean(rel_diff_single)) %>% 
  mutate(parameter_label = case_when( ## recode parameter values
    parameter == "rho" ~ "RBC replenishment (ρ)",
    parameter == "phin" ~ "Half-life indis (ϕn)",
    parameter == "phiw" ~ "Half-life targeted (ϕw)",
    parameter == "psin" ~ "Activation indis (ψn)",
    parameter == "psiw" ~ "Activation targeted (ψw)",
    parameter == "beta" ~ "Burst size (β)",
  ))

`summarise()` has grouped output by 'long_label'. You can override using the `.groups` argument.

## for violin plots, what we can do is plot out only the 89% credible interval so the graph is easier to interpret
mc_single_fitness.long_p_f <- mc_single_fitness.long_p %>% 
  left_join(select(mc_single_fitness.sum, long_label, parameter, ci_lower, ci_higher, parameter_label), by = c("long_label", "parameter")) %>% 
  filter(rel_diff_single >= ci_lower & rel_diff_single < ci_higher)


## arrange ordering of cues and parameters
mc_single_fitness.long_p_f$long_label <- factor(mc_single_fitness.long_p_f$long_label,
                                           c("RBC log &\nasexual iRBC log",
                                             "Asexual iRBC", "Asexual iRBC log",
                                             "Sexual iRBC", "Sexual iRBC log",
                                             "Total iRBC", "Total iRBC log",
                                             "Gametocyte", "Gametocyte log", 
                                             "RBC", "RBC log"))

mc_single_fitness.long_p_f$parameter_label <- factor(mc_single_fitness.long_p_f$parameter_label,
                                                     c("Burst size (β)",
                                                       "RBC replenishment (ρ)",
                                                       "Half-life indis (ϕn)",
                                                       "Half-life targeted (ϕw)",
                                                       "Activation indis (ψn)",
                                                       "Activation targeted (ψw)"))

plotting 89% credible interval. By visual inspection of the distribution, ci represented the distribution (more honestly) than quantile, even when quantile seems to display larger differences of logging cues.

#——— arrange plots together ———#

ggarrange(mc_all_fitness.pl, mc_partition.pl, align = "h", widths = c(1.3, 2), labels = c("A", "B"))
ggsave(units = "px", dpi = 300, width = 2250, height = 1500, filename = here("code_repository/figures/mc_fitness_partition.tiff"), bg = "white", scale = 1.2)

#==========================================# # Supplementary figure on MC posterior parameter distribution #==========================================# ## Process data

## manually make dataframe containing the parameter values used in the deterministic model
par_det.df <- data.frame(
  parameter = c("rho", "phi_N1", "phi_N2", "iota_N1", "iota_N2", "burst"),
  deterministic = c(2.627156e-01, 3.520591e-02, 5.508420e+02, 1.669234e+01, 8.431785e-01, 5.721000e+00 )
)

## make into long format
posterior.long <- posterior.df %>% 
  tidyr::pivot_longer(-id, names_to = "parameter") %>% 
  left_join(par_det.df, by = "parameter") %>% ## add in the deterministic parameter values
  mutate(label = case_when(
    parameter == "rho" ~ "RBC replenishment (ρ)",
    parameter == "phi_N1" ~ "Half-life indis (ϕn)",
    parameter == "phi_N2" ~ "Half-life targeted (ϕw)",
    parameter == "iota_N1" ~ "Activation indis (ψn)",
    parameter == "iota_N2" ~ "Activation targeted (ψw)",
    parameter == "burst" ~ "Burst size (β)"
  )) ## rename paramter values

plot

## portion of parameter values that does not need log-transforming
posterior_1.pl <- ggplot(posterior.long %>% filter(!parameter %in% c("phi_N1", "phi_N2", "iota_N1"))) +
  geom_density(aes(x = value), fill = "grey") +
  geom_vline(aes(xintercept = deterministic), linetype = "dashed") +
  facet_wrap(~label, scales = "free") +
  labs(x = "", y = "Density") +
  theme_bw() +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank())

## Plots where the x-axis should be log-transformed so it is easier to read
posterior_2.pl <- ggplot(posterior.long %>% filter(parameter %in% c("phi_N1", "phi_N2", "iota_N1"))) +
  geom_density(aes(x = value), fill = "grey") +
  geom_vline(aes(xintercept = deterministic), linetype = "dashed") +
  facet_wrap(~label, scales = "free") +
  labs(x = "Value", y = "Density") +
  scale_x_continuous(trans = "log10") +
  theme_bw() +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank())

## plot together
ggarrange(posterior_1.pl, posterior_2.pl, ncol = 1, align = "hv")

Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y,  :
  conversion failure on 'Activation targeted (ψw)' in 'mbcsToSbcs': dot substituted for <cf>
Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y,  :
  conversion failure on 'Activation targeted (ψw)' in 'mbcsToSbcs': dot substituted for <88>
Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y,  :
  conversion failure on 'Burst size (β)' in 'mbcsToSbcs': dot substituted for <ce>
Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y,  :
  conversion failure on 'Burst size (β)' in 'mbcsToSbcs': dot substituted for <b2>
Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y,  :
  conversion failure on 'RBC replenishment (ρ)' in 'mbcsToSbcs': dot substituted for <cf>
Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y,  :
  conversion failure on 'RBC replenishment (ρ)' in 'mbcsToSbcs': dot substituted for <81>
Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y,  :
  conversion failure on 'Activation targeted (ψw)' in 'mbcsToSbcs': dot substituted for <cf>
Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y,  :
  conversion failure on 'Activation targeted (ψw)' in 'mbcsToSbcs': dot substituted for <88>
Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y,  :
  conversion failure on 'Burst size (β)' in 'mbcsToSbcs': dot substituted for <ce>
Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y,  :
  conversion failure on 'Burst size (β)' in 'mbcsToSbcs': dot substituted for <b2>
Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y,  :
  conversion failure on 'RBC replenishment (ρ)' in 'mbcsToSbcs': dot substituted for <cf>
Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y,  :
  conversion failure on 'RBC replenishment (ρ)' in 'mbcsToSbcs': dot substituted for <81>
Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y,  :
  conversion failure on 'Activation targeted (ψw)' in 'mbcsToSbcs': dot substituted for <cf>
Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y,  :
  conversion failure on 'Activation targeted (ψw)' in 'mbcsToSbcs': dot substituted for <88>
Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y,  :
  conversion failure on 'Burst size (β)' in 'mbcsToSbcs': dot substituted for <ce>
Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y,  :
  conversion failure on 'Burst size (β)' in 'mbcsToSbcs': dot substituted for <b2>
Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y,  :
  conversion failure on 'RBC replenishment (ρ)' in 'mbcsToSbcs': dot substituted for <cf>
Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y,  :
  conversion failure on 'RBC replenishment (ρ)' in 'mbcsToSbcs': dot substituted for <81>
Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y,  :
  conversion failure on 'Activation targeted (ψw)' in 'mbcsToSbcs': dot substituted for <cf>
Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y,  :
  conversion failure on 'Activation targeted (ψw)' in 'mbcsToSbcs': dot substituted for <88>
Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y,  :
  conversion failure on 'Burst size (β)' in 'mbcsToSbcs': dot substituted for <ce>
Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y,  :
  conversion failure on 'Burst size (β)' in 'mbcsToSbcs': dot substituted for <b2>
Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y,  :
  conversion failure on 'RBC replenishment (ρ)' in 'mbcsToSbcs': dot substituted for <cf>
Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y,  :
  conversion failure on 'RBC replenishment (ρ)' in 'mbcsToSbcs': dot substituted for <81>
Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y,  :
  conversion failure on 'Activation indis (ψn)' in 'mbcsToSbcs': dot substituted for <cf>
Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y,  :
  conversion failure on 'Activation indis (ψn)' in 'mbcsToSbcs': dot substituted for <88>
Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y,  :
  conversion failure on 'Half-life indis (ϕn)' in 'mbcsToSbcs': dot substituted for <cf>
Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y,  :
  conversion failure on 'Half-life indis (ϕn)' in 'mbcsToSbcs': dot substituted for <95>
Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y,  :
  conversion failure on 'Half-life targeted (ϕw)' in 'mbcsToSbcs': dot substituted for <cf>
Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y,  :
  conversion failure on 'Half-life targeted (ϕw)' in 'mbcsToSbcs': dot substituted for <95>
Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y,  :
  conversion failure on 'Activation indis (ψn)' in 'mbcsToSbcs': dot substituted for <cf>
Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y,  :
  conversion failure on 'Activation indis (ψn)' in 'mbcsToSbcs': dot substituted for <88>
Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y,  :
  conversion failure on 'Half-life indis (ϕn)' in 'mbcsToSbcs': dot substituted for <cf>
Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y,  :
  conversion failure on 'Half-life indis (ϕn)' in 'mbcsToSbcs': dot substituted for <95>
Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y,  :
  conversion failure on 'Half-life targeted (ϕw)' in 'mbcsToSbcs': dot substituted for <cf>
Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y,  :
  conversion failure on 'Half-life targeted (ϕw)' in 'mbcsToSbcs': dot substituted for <95>
Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y,  :
  conversion failure on 'Activation indis (ψn)' in 'mbcsToSbcs': dot substituted for <cf>
Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y,  :
  conversion failure on 'Activation indis (ψn)' in 'mbcsToSbcs': dot substituted for <88>
Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y,  :
  conversion failure on 'Half-life indis (ϕn)' in 'mbcsToSbcs': dot substituted for <cf>
Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y,  :
  conversion failure on 'Half-life indis (ϕn)' in 'mbcsToSbcs': dot substituted for <95>
Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y,  :
  conversion failure on 'Half-life targeted (ϕw)' in 'mbcsToSbcs': dot substituted for <cf>
Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y,  :
  conversion failure on 'Half-life targeted (ϕw)' in 'mbcsToSbcs': dot substituted for <95>
Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y,  :
  conversion failure on 'Activation indis (ψn)' in 'mbcsToSbcs': dot substituted for <cf>
Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y,  :
  conversion failure on 'Activation indis (ψn)' in 'mbcsToSbcs': dot substituted for <88>
Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y,  :
  conversion failure on 'Half-life indis (ϕn)' in 'mbcsToSbcs': dot substituted for <cf>
Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y,  :
  conversion failure on 'Half-life indis (ϕn)' in 'mbcsToSbcs': dot substituted for <95>
Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y,  :
  conversion failure on 'Half-life targeted (ϕw)' in 'mbcsToSbcs': dot substituted for <cf>
Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y,  :
  conversion failure on 'Half-life targeted (ϕw)' in 'mbcsToSbcs': dot substituted for <95>

  
## save
ggsave(units = "px", dpi = 300, width = 2000, height = 1300, filename = here("code_repository/figures/posterior.tiff"), bg = "white", scale = 1)

#==========================================# # Rank plot of parasite fitness with only one parameter varying #==========================================# ## Process data

## get median fitness values
mc_single_fitness.sum2 <- mc_single_fitness.long2 %>% 
  group_by(short_label, parameter_label) %>% 
  summarize(median = median(value),
            mean = mean(value),
            geo_mean = exp(mean(log(value))))

`summarise()` has grouped output by 'short_label'. You can override using the `.groups` argument.

plot

#==================================# # Deciphering why certain cues are less susceptible to burst size variation #==================================# #————— Optimized vs non-optimized fitness against burst size variation————-# This could be attributed to lower ability of cues (even with optimal reaction norm) to respond to changed beta. Or it could be robustness! ## Process data for downstream analysis (execute once)

## get list of files
mc_burst.ls <- list.files(here("code_repository/data/mc_burst_opt/"), pattern = "*.csv", full.names = T)

## read in and bind. Note that we are excluding the parameter values for now because R log + I log has 9 of them
mc_burst.df_raw <- do.call(rbind, lapply(mc_burst.ls, function(x){
  df <- read.csv(x)
  df_p <- df %>% select(id, cue, log, beta, fitness, default)
  return(df_p)}))

## we note that the R log + I log has a slightly higher fitness than time (df = 3), which is due to the greater flexibility of spline. NOte in previous optimization, increasing the df = 9 allowed time based cue to achieve higher fitness than R log + I log
mc_burst.df_raw %>% filter(id == "R_log+I_log")

## we optimized each cue (and burst size) via 2 methods: one starting from the optimal parameter size when beta = 5.721 (default = F) and one starting at (0.5x4) (default = T). We are going to pick the whichever one of them gave the highest fitness
mc_burst.df <- mc_burst.df_raw %>% 
  filter(id != "time") %>% ## time is ignored here because irrelavance
  group_by(id, beta) %>% 
  top_n(1, fitness) %>% ### pick the one with the highest fitness
  mutate(log = ifelse(log == "log10", "log", "none")) ## change log classification for consistency

## write 
write.csv(mc_burst.df, here("code_repository/data/mc_burst_opt.csv"))

Process data for plotting

## list of burst size used in the optimization simulations
burst_ls <- c(3.98,4.7,5.53, 6.08)

## here, we want to see if after optimization, the optimized fitness values correlated with unoptimized fitness values. That is, is the decline of fitness caused by burst size variation due to limitations of the cue to achieve higher fitness or robustness of the cue?
posterior_burst_iter <- posterior.df %>% 
  filter(round(burst, digits = 2) %in% burst_ls) %>% ## get iterations that has a burst value close to the one we tested
  mutate(diff_3 = abs(burst - burst_ls[[1]]),
         diff_4 = abs(burst - burst_ls[[2]]),
         diff_5 = abs(burst - burst_ls[[3]]),
         diff_6 = abs(burst - burst_ls[[4]]))

## get iterations that we should take fitness values from
burst_iter_3 <- posterior_burst_iter[posterior_burst_iter$diff_3 == min(posterior_burst_iter$diff_3), "id"]
burst_iter_4 <- posterior_burst_iter[posterior_burst_iter$diff_4 == min(posterior_burst_iter$diff_4), "id"]
burst_iter_5 <- posterior_burst_iter[posterior_burst_iter$diff_5 == min(posterior_burst_iter$diff_5), "id"]
burst_iter_6 <- posterior_burst_iter[posterior_burst_iter$diff_6 == min(posterior_burst_iter$diff_6), "id"]

burst_iter_ls <- c(burst_iter_3, burst_iter_4, burst_iter_5, burst_iter_6)

## get joining dataframe to associate iteration number of burst size
burst_iter.df <- posterior.df %>% 
  filter(id %in% c( burst_iter_ls)) %>% 
  mutate(burst_round = round(burst, digits = 2)) %>% 
  select(iter = id, beta = burst_round)

mc_burst_fitness_final.df <- mc_single_fitness.long %>% 
  filter(parameter == "beta" & iter %in% c(burst_iter_ls)) %>% ## keep only rows where the parameter value altered is beta and belongs to the appropriate iterations 
  left_join(burst_iter.df, by = "iter") %>% ## get beta value associated with iterations
  select(cue, log, fitness_unopt = value, beta) %>% ## select relevant columns
  left_join(select(mc_burst.df, cue, log, beta, fitness_opt = fitness), by = c("cue", "log", "beta")) %>%  ## join to get fitness value after optiization based on cue, log status, and beta value
  left_join(select(ez_label, cue, log, short_label, long_label), by = c("cue", "log")) %>% 
  mutate(short_label = ifelse(cue == "R log & I log", cue, short_label)) %>%  ## manually generate label for R log and I log
  pivot_longer(cols = c(fitness_unopt, fitness_opt)) %>% ## make long
  #select(-id) %>% 
  mutate(classification = ifelse(name == "fitness_unopt", "Unoptimized", "Optimized"))

## get mean values
mc_burst_fitness_med.df <- mc_burst_fitness_final.df %>% 
  group_by(classification, short_label) %>% 
  mutate(mean = mean(value),
         geom_mean = exp(mean(log(value)))) %>%
  distinct(cue, log, .keep_all = T) %>% 
  ungroup()

plot

## note that all points are ranked by increasing geometric mean within each facets!
mc_burst_fitness.pl <- ggplot() +
  geom_line(data = mc_burst_fitness_final.df, aes(x = tidytext::reorder_within(short_label, value, classification, fun = function(x) exp(mean(log(x)))), 
                                                  y = value, group = beta), color = "dark grey") +
  geom_point(data = mc_burst_fitness_final.df, aes(x = tidytext::reorder_within(short_label, value, classification, fun = function(x) exp(mean(log(x)))), 
                                                   y = value, color = beta), size = 2.5, alpha = 0.6) +
  geom_line(data = mc_burst_fitness_med.df, aes(x = tidytext::reorder_within(short_label, value, classification, fun = function(x) exp(mean(log(x)))), 
                                                y = geom_mean, group = classification), size = 1.5) +
  geom_point(data = mc_burst_fitness_med.df, aes(x = tidytext::reorder_within(short_label, value, classification, fun = function(x) exp(mean(log(x)))), 
                                                 y = geom_mean), size = 3, shape = 17, color = "black") +
  facet_wrap(~classification, ncol = 1, scales = "free") +
  tidytext::scale_x_reordered() +
  scale_color_viridis_c() +
  labs(x = "Cue(s)", y = "Fitness", color = "Burst size") +
  theme_bw() +
  theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.text.x = element_text(angle = 45, hjust=1))

#————— The advantage of I log is not due to default optimization position—————-# Here, we will use the 4 optimized parameter at various time point and obtain their fitness across the 3 other burst values ## function to get fitness and conversion rate

get_mc_fitness_diff <- function(df, dual = F){
  ## list of burst values to test
  burst_ls <- c(3.98,4.7,5.53, 5.721, 6.08)
  
  ## get parameter set
  ### assign the baseline parameter set
  parameters_tsukushi_default <- c(R1 = 8.89*10^6, 
                lambda = 3.7*10^5,
                mu = 0.025, 
                p = 8*10^-6, 
                alpha = 1, 
                alphag = 2, 
                beta = df$beta, ## beta is inserted here (place holder)
                mum = 48, 
                mug = 4, 
                I0 = 43.85965, 
                Ig0 = 0, 
                a = 150, 
                b = 100, 
                sp = 1,
                psin = 16.69234,
                psiw = 0.8431785,
                phin = 0.03520591, 
                phiw = 550.842,
                iota = 2.18*(10^6),
                rho = 0.2627156)
  
  ### get list of paramaters with all 4 values of beta
  parameter_tsukushi_ls <- lapply(burst_ls, function(x) replace(parameters_tsukushi_default, 7, x))
  
  ## read in input for model simulation
  if(isFALSE(dual)){
    par <- c(df$par1, df$par2, df$par3, df$par4) ## parameter set
    cue <- df$cue
    log <- df$log
    cue_range <- seq(df$low, df$high, df$by)
  } else{
    par <- c(df$par1, df$par2, df$par3, df$par4, df$par5, df$par6, df$par7, df$par8, df$par9)
    cue <- "R"
    cue_b <- "I"
    log <- "log10"
    log_b <- "log10"
    cue_range <- seq(6, 7, length.out = 500)
    cue_range_b <- seq(0,   6.77815125, length.out = 500)
  }
  
  ## get dynamics data
  if(isFALSE(dual)){
    fitness_ls <- lapply(parameter_tsukushi_ls, function(x){
      dyn <- chabaudi_si_clean(parameters_cr = par, 
      parameters = x, 
      immunity = "tsukushi",
      time_range = seq(0, 20, 0.01), 
      cue = cue, 
      log_cue = log, 
      cue_range = cue_range, 
      solver = "vode",
      dyn = T)
      ### process dynamics dataframe to add beta, cue, and log
      dyn_p <- dyn %>% mutate(id = df$id, beta_optimized = df$beta, beta_simulated = x[[7]])
      write_parquet(dyn_p, here(paste0("code_repository/data/mc_burst_dyn/", df$id, "_", gsub("\\.", "-", df$beta),
                                       "_", gsub("\\.", "-", x[[7]]), "_dyn.parquet")))
      ### get fitness
      fitness <- max(dyn[dyn$variable == "tau_cum", 3], na.rm = T)
      fitness_beta <- c(beta_simulated = x[[7]], fitness = fitness) ## attach the beta used in the simulation here
      return(fitness_beta)
    })
  } else{ ## dual cue models
    fitness_ls <- lapply(parameter_tsukushi_ls, function(x){
      dyn <- chabaudi_si_clean(parameters_cr = par, 
      parameters = x,
      immunity = "tsukushi",
      time_range = seq(0, 20, 0.01), 
      cue = cue, 
      cue_b = cue_b,
      log_cue = log,
      log_cue_b = log_b,
      cue_range = cue_range, 
      cue_range_b = cue_range_b,
      gam = "te",
      solver = "vode",
      dyn = T)
      
     ### process dynamics dataframe to add beta, cue, and log
      dyn_p <- dyn %>% mutate(id = df$id, beta_optimized = df$beta, beta_simulated = x[[7]])
      write_parquet(dyn_p, here(paste0("code_repository/data/mc_burst_dyn/", df$id, "_", gsub("\\.", "-", df$beta),
                                       "_", gsub("\\.", "-", x[[7]]), "_dyn.parquet")))
      ### get fitness
      fitness <- max(dyn[dyn$variable == "tau_cum", 3], na.rm = T)
      
      fitness_beta <- c(beta_simulated = x[[7]], fitness = fitness) ## attach the beta used in the simulation here
      return(fitness_beta)
      }
    )
  }
  
  ## add identifiers to results (fitness)
  res_final <- do.call(rbind, fitness_ls) %>% 
    as.data.frame() %>% 
    mutate(id = df$id,
           cue = cue,
           log = log,
           beta_optimized = df$beta)
  
  return(res_final)
}

get input dataframe for function to run

## get list of files
mc_burst.ls <- list.files(here("code_repository/data/mc_burst_opt/"), pattern = "*.csv", full.names = T)
mc_RI_burst.ls <- list.files(here("code_repository/data/mc_burst_opt/"), pattern = "R_log_I_log*", full.names = T)

## read in and bind. Exclude R log + I log due to differences in # of parameters
mc_burst_par.df_raw <- do.call(rbind, lapply(mc_burst.ls[!grepl("R_log_I_log", mc_burst.ls)], function(x){
  df <- read.csv(x)
  df_p <- df %>% select(id, cue, log, beta, fitness, par1, par2, par3, par4, default)
  return(df_p)}))
mc_burst_RI.df_raw <- do.call(rbind, lapply(mc_RI_burst.ls, function(x){
  df <- read.csv(x)
  df_p <- df %>% select(id, cue, log, beta, fitness, par1, par2, par3, par4, par5, par6, par7, par8, par9, default)
  return(df_p)}))

## get the top strategy based on fitness
mc_burst_par.df_f <- mc_burst_par.df_raw %>% 
  group_by(id, beta) %>% 
  top_n(1, fitness) %>% 
  filter(id != "time") ## get rid of time

mc_burst_RI.df <- mc_burst_RI.df_raw %>% 
  group_by(id, beta) %>% 
  top_n(1, fitness)

## join with si_opt to get cue ranges
mc_burst_par.df_fp <- mc_burst_par.df_f %>% 
  left_join(select(si_opt.df, id, low, high, by, long_label), by = "id")

## write
write.csv(mc_burst_par.df_fp, here("code_repository/data/mc_burst_single_input.csv"))
write.csv(mc_burst_RI.df, here("code_repository/data/mc_burst_dual_input.csv"))

run function

process data for plotting

## get geometric mean. note that we are including cases where optimized beta = simulated beta to account for how good the rn is under optimal circumstances
mc_burst_diff.sum <- mc_burst_final_diff.df %>% 
  group_by(short_label) %>% 
  summarise(geom_mean = exp(mean(log(fitness)))) %>% 
  arrange(geom_mean)

## reorder x axis based on geometric mean
mc_burst_final_diff.df$short_label <- factor(mc_burst_final_diff.df$short_label, mc_burst_diff.sum$short_label)

plot

mc_burst_diff.pl <- ggplot() +
  geom_violin(data = mc_burst_final_diff.df, 
             aes(x = short_label, y = fitness), color = "transparent", fill = "light grey") +
  geom_jitter(data = mc_burst_final_diff.df, 
             aes(x = short_label, y = fitness, color = beta_optimized), size = 2, alpha = 0.7, width = 0.2) +
  scale_color_viridis_c() +
  geom_crossbar(data = mc_burst_diff.sum, aes(
    x = short_label, y = geom_mean,
    xmin = short_label, xmax = short_label, ymin = geom_mean, ymax = geom_mean),
                  size=1,col="black", width = .5) +
  labs(x = "Cue(s)", y = "Fitness", color = "Optimized\nburst size") +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 45, hjust=1))

#————— Changes in reaction norm in response to burst size variation————-# ## Function to obtain dynamics data (rug) for parasites adopting the optimal strategy for each burst size (execute once)

get_mc_burst_rug <- function(df, dual = F){
  ## get paramater set
  parameters_tsukushi_burst <- c(R1 = 8.89*10^6, 
                lambda = 3.7*10^5,
                mu = 0.025, 
                p = 8*10^-6, 
                alpha = 1, 
                alphag = 2, 
                beta = df$beta, ## beta is inserted here
                mum = 48, 
                mug = 4, 
                I0 = 43.85965, 
                Ig0 = 0, 
                a = 150, 
                b = 100, 
                sp = 1,
                psin = 16.69234,
                psiw = 0.8431785,
                phin = 0.03520591, 
                phiw = 550.842,
                iota = 2.18*(10^6),
                rho = 0.2627156)
  
  ## read in input
  if(isFALSE(dual)){
    par <- c(df$par1, df$par2, df$par3, df$par4) ## parameter set
    cue <- df$cue
    log <- df$log
    cue_range <- seq(df$low, df$high, df$by)
  } else{
    par <- c(df$par1, df$par2, df$par3, df$par4, df$par5, df$par6, df$par7, df$par8, df$par9)
    cue <- "R"
    cue_b <- "I"
    log <- "log10"
    log_b <- "log10"
    cue_range <- seq(6, 7, length.out = 500)
    cue_range_b <- seq(0,   6.77815125, length.out = 500)
  }
  
  ## get dynamics data
  if(isFALSE(dual)){
    dyn <- chabaudi_si_clean(parameters_cr = par, 
    parameters = parameters_tsukushi_burst, 
    immunity = "tsukushi",
    time_range = seq(0, 20, 0.01), 
    cue = cue, 
    log_cue = log, 
    cue_range = cue_range, 
    dyn = T)
  } else{
    dyn <- chabaudi_si_clean(parameters_cr = par, 
    parameters = parameters_tsukushi_burst, 
    immunity = "tsukushi",
    time_range = seq(0, 20, 0.01), 
    cue = cue, 
    cue_b = cue_b,
    log_cue = log,
    log_cue_b = log_b,
    cue_range = cue_range, 
    cue_range_b = cue_range_b,
    gam = "te",
    dyn = T)
  }
  
  ## get rug data only. Note that for dual cue models we are filtering out R and I 
  ### filter to get only the relevant columns
  if(isFALSE(dual)){
    if(cue != "I+Ig"){
    dyn_f <- dyn %>% filter(variable == cue)
  } else{
    dyn_f <- dyn %>% 
      filter(variable %in% c("I", "Ig")) %>% 
      group_by(time) %>% 
      summarise(value = sum(value, na.rm = T)) %>% 
      mutate(variable = "I+Ig")
    }
  } else{
    dyn_f <- dyn %>% filter(variable %in% c("R", "I"))
  }
  
  ### log-transform if needed. This does not needed to be ammended for dual cue models because both value is logged in the R log & I log model. Please change this if working with incidences where only one of the variable is log-transformed
  if(log == "log10"){
    dyn_fp <- dyn_f %>% 
      mutate(value = log10(value),
             value = ifelse(value == -Inf, 0, value)) ## for -Inf (which is what happens when you log transform 0), reset to 0
  } else{
    dyn_fp <- dyn_f
  }
  
  ## add identifiers to results
  dyn_final <- dyn_fp %>% 
    mutate(id = df$id,
           cue = cue,
           log = log,
           beta = df$beta)
  
  return(dyn_final)
}

Prepare rug and reaction norm (execute once)

mc_burst_rn_RI.df_p <- mc_burst_rn_RI.df %>% 
  left_join(mc_burst_bound_RI %>% filter(variable == "R"), by = c("id", "beta")) %>% 
  left_join(mc_burst_bound_RI %>% filter(variable == "I") %>% select(id, beta, min_b = min, max_b = max), by = c("id", "beta")) %>% 
  filter(cue_range >= min & cue_range <= max & cue_range_b >= min_b & cue_range_b <= max_b)

Adding missing grouping variables: `variable`

Prepare data for plotting

## get label and only look at I log10 and I log + R log10 for comparison
mc_burst_dyn.df_p <- mc_burst_dyn.df %>% 
  filter(id %in% c("I_log","R_log+I_log") & variable == "cr" & beta_simulated != 5.721) %>% ## only get the cues we want to plot
  left_join(ez_label, by = "id") %>% 
  mutate(long_short_label = ifelse(id == "R_log+I_log", "RBC log & asexual iRBC log\n(R log & I log)",
    paste0(long_label, " (", short_label, ")")))

Plot

ggplot() +
  geom_line(data = mc_burst_dyn.df_p %>% filter(beta_optimized != beta_simulated), 
            aes(x = time, y = value, group = beta_optimized, color = beta_optimized)) + ## only plotting unoptimized
  geom_line(data = mc_burst_dyn.df_p %>% filter(beta_optimized == beta_simulated), 
            aes(x = time, y = value, group = beta_optimized, color = beta_optimized), linetype = "longdash", size = 1) + ## optimized cr is indicated by a dashed line
  facet_grid(rows = vars(beta_simulated), cols = vars(long_short_label)) +
  labs(x = "Time (days)", y = "Conversion rate", color = "Optimized\nburst size") +
  scale_color_viridis_c() +
  theme_bw() +
  theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank())

Warning: Removed 4 rows containing missing values (`geom_line()`).
Warning: Removed 4 rows containing missing values (`geom_line()`).

#—— arrange together!————-#

## first arrange the two fitness graphs
mc_burst_fitness_com.pl <- ggarrange(mc_burst_fitness.pl, mc_burst_diff.pl, ncol = 1, align = "hv", labels = c("A", "B"), heights = c(0.6,0.4))

Warning: Graphs cannot be horizontally aligned unless the axis parameter is set. Placing graphs unaligned.

## put them together
ggarrange(mc_burst_fitness_com.pl, mc_burst_cr.pl, ncol = 2, align = "hv", labels = c("", "C"), widths = c(0.5,0.55))

Warning: Removed 4 rows containing missing values (`geom_line()`).
Warning: Removed 4 rows containing missing values (`geom_line()`).
Warning: Graphs cannot be vertically aligned unless the axis parameter is set. Placing graphs unaligned.
Warning: Graphs cannot be horizontally aligned unless the axis parameter is set. Placing graphs unaligned.

ggsave(units = "px", dpi = 300, width = 2650, height = 1500, filename = here("code_repository/figures/mc_burst_main.tiff"), bg = "white", scale = 1.26)

#=================================# # Supplementary fig: all rn at different burst sizes #=================================# ## Process data for plotting

## get plotting label
mc_burst_rn.df_p <- mc_burst_rn.df %>% 
  left_join(select(ez_label, id, short_label), by = "id") %>% 
  mutate(long_short_label = paste0(long_label, " (", short_label, ")"))

## get plotting order. We are ranking them based on the geometric mean of unoptimized fitness (for all combinations)
mc_burst_order <- mc_burst_final_diff.df %>% 
  mutate(long_label = ifelse(id == "R_log+I_log", "RBC log & asexual iRBC log", long_label)) %>% 
  mutate(long_short_label = paste0(long_label, " (", short_label, ")")) %>% 
  group_by(long_short_label) %>% 
  summarise(geom_fitness = exp(mean(log(fitness)))) %>% 
  ungroup() %>%
  arrange(-geom_fitness)

mc_burst_rn.df_p$long_short_label <- factor(mc_burst_rn.df_p$long_short_label, mc_burst_order$long_short_label)

## Normalize cue range for plotting
mc_burst_rn.df_p_n <- mc_burst_rn.df_p %>% 
  group_by(id) %>% 
  mutate(range = max(cue_range) - min(cue_range),
         mean = mean(cue_range),
         cue_range_norm = (cue_range - mean)/range)

## make rug plot of RI wider for trajectory plotting
mc_burst_RI_rug.wide <- mc_burst_RI_rug.df %>% 
  select(time, variable,value, beta) %>% 
  pivot_wider(names_from = variable, values_from = value, id_cols = c(time, beta))

mc_burst_RI_rug.wide

plot

## single cue
mc_burst_rn_sc.pl <- ggplot(data = mc_burst_rn.df_p_n) +
  geom_line(aes(x = cue_range_norm, y = cr, color = beta, group = beta)) +
  facet_wrap(~long_short_label, scales = "free", ncol = 3) +
  labs(x = "Cue range (normalized)", y = "Conversion rate", color = "Burst size") +
  scale_color_viridis_c() +
  theme_bw() +
  theme(legend.position = "bottom",
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank()) +
  scale_x_continuous(labels = label_number(accuracy = 0.1))

## dual cue
mc_burst_rn_dc.pl <- ggplot() +
  geom_raster(data = mc_burst_rn_RI.df, aes(y = cue_range, x = cue_range_b, fill = cr)) +
  geom_path(data = mc_burst_RI_rug.wide, aes(x = I, y = R), size = 0.75, color = "cyan", arrow = arrow(angle = 30, length = unit(0.1, "inches"))) +
  geom_point(data = mc_burst_RI_rug.wide %>% filter(time %% 1 == 0), aes(x = I, y = R), size = 1, color = "cyan") +
  scale_fill_viridis_c(option = "C") +
  facet_wrap(~beta, ncol = 1) +
  labs(x = "Asexual iRBC log10", y = "RBC log10", fill = "Conversion\nrate") +
  theme_bw() +
  theme(legend.position = "bottom",
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
    plot.subtitle=element_text(size=13)) +
  scale_y_continuous(labels = label_number(accuracy = 0.1),
                     limits = c(6.8,7))

## plot together
ggarrange(mc_burst_rn_sc.pl, mc_burst_rn_dc.pl, widths = c(2.5,1.1), labels = c("A", "B"), ncol = 2, nrow = 1, align = "hv")
ggsave(units = "px", dpi = 300, width = 2550, height = 1800, filename = here("code_repository/figures/mc_burst_rn.tiff"), bg = "white", scale = 1)

sanity check. RI rn is different between the beta != 3.98

mc_burst_rn_RI.df %>% 
  filter(beta != 3.98) %>% 
  group_by(cue_range, cue_range_b) %>% 
  mutate(sd = sd(cr)) %>% 
  ggplot() + geom_raster(aes(x = cue_range_b, y = cue_range, fill = sd))

#=====================================# # Impact of cue perception on virulence #====================================# #——— Divid dynamics data into high and low performing —————# # Single cue

top_4_single_cue$short_label

[1] "Ig log"   "I+Ig log" "G log"    "I log"   

Dual cue

dual_cue_dyn_30.df %>% group_by(id, id_b) %>% summarise(max_time = max(time))

`summarise()` has grouped output by 'id'. You can override using the `.groups` argument.

#————-Plot disease curves—————–# # Single cue

dual cue

#———area of curve vs fitness—————# # function to calculate area within curve

single cue area

dual cue

fit GLM

# fit lm
si_lm <- lm(fitness_norm ~ area, data = si_opt_area.df)
dual_lm <- lm(fitness_norm ~ area, data = dual_cue_f_final_area.df)

# get prediction with 95% CI
si_lm.ci = predict(si_lm, interval = "confidence") %>% cbind.data.frame(area = si_opt_area.df$area)
dual_lm.ci = predict(dual_lm, interval = "confidence") %>% cbind.data.frame(area = si_opt_area.df$area)

# look at R^2
summary(si_lm)$r.squared ## 0.1765248

[1] 0.1765248

summary(dual_lm)$r.squared ## 0.1053591

[1] 0.1053591

plot correlation

# single cue
si_corr_area.plt <- ggplot() +
  geom_point(data = si_opt_area.df, aes(x = area, y = fitness_norm, color = "Single cue", shape = "Single cue"), size = 3, alpha = 0.7) +
  geom_line(data = si_lm.ci, aes(x = area, fit)) +
  geom_ribbon(data = si_lm.ci, aes(x = area, ymin = lwr, ymax = upr), alpha = .15) +
  annotate(geom = "text", label = "R^2=0.18", x = 1.6*(10^11), y = 0.76) +
  labs(x = "Area", y = "Normalized fitness", shape = "Cue category", color = "Cue category") +
  scale_color_manual(values = c(orange)) +
  xlim(min(si_opt_area.df$area), max(dual_cue_f_final_area.df$area)) +
  theme_classic() +
  theme(legend.position = "none")

# dual cue
dual_corr_area.plt <- ggplot() +
  geom_point(data = dual_cue_f_final_area.df, aes(x = area, y = fitness_norm, color = "Dual cue", shape = "Dual cue"), size = 3, alpha = 0.7) +
  geom_line(data = dual_lm.ci, aes(x = area, fit)) +
  geom_ribbon(data = dual_lm.ci, aes(x = area, ymin = lwr, ymax = upr), alpha = .15) +
 annotate(geom = "text", label = "R^2=0.11", x = 1.6*(10^11), y = 0.76) +
  labs(x = "Area", y = "Normalized fitness", shape = "Cue category", color = "Cue category") +
  scale_color_manual(values = c(blue)) +
  xlim(min(si_opt_area.df$area), max(dual_cue_f_final_area.df$area)) +
  theme_classic() +
  theme(legend.position = "none")

#———assemble figure————#

si_dc_area.plt <- ggarrange(si_dc.plt, si_corr_area.plt, ncol = 2, nrow = 1, align = "hv", common.legend = T, labels = c("A", "B"))
dual_dc_area.plt <- ggarrange(dual_dc.plt, dual_corr_area.plt, ncol = 2, nrow = 1, align = "hv", common.legend = T, labels = c("C", "D"))

Warning: Removed 2 rows containing missing values or values outside the scale range (`geom_point()`).
Warning: Removed 2 rows containing missing values or values outside the scale range (`geom_line()`).

ggarrange(si_dc_area.plt, dual_dc_area.plt, nrow = 2, align = "v")

ggsave(units = "px", dpi = 300, width = 2550, height = 2000, filename = here("code_repository/figures/virulence.tiff"), bg = "white", scale = 1)

#=============================# # Validation of optimization #=============================# #—— process input data ——# # Process data

## get difference betweent fitness conferred by model produced by a random spline strategy "V1" vs optimized fitness
validation.df_p <- validation.df %>% 
  left_join(select(si_opt.df, id, fitness_20, short_label, long_label), by = c("id")) %>% 
  mutate(diff = fitness_20-V1,
         long_short_label = paste0(long_label, " (", short_label, ")"))

## no random strategy performed better!
validation.df_p %>% filter(diff<=0)

Plot

ggplot(validation.df_p) +
  geom_density(aes(x = diff), fill = "light grey") +
  facet_wrap(~long_short_label, scales = "free") +
  labs(x = "Fitness difference (Optimal-Random)", y = "Density") +
  theme_bw() +
  theme(legend.position = "bottom",
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank()) 

ggsave(units = "px", dpi = 300, width = 2550, height = 1500, filename = here("code_repository/figures/si_validation.tiff"), bg = "white", scale = 1)